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Standard Ewald sums, which calculate e.g. the electrostatic energy or the force in periodically closed 
systems of charged particles, can be efficiently speeded up by the use of the Fast Fourier Transfor- 
mation (FFT). In this article we investigate three algorithms for the FFT-accelerated Ewald sum, 
which attracted a widespread attention, namely, the so-called particle-particle-particle-mesh (P 3 M), 
particle mesh Ewald (PME) and smooth PME method. We present a unified view of the underlying 
techniques and the various ingredients which comprise those routines. Additionally, we offer detailed 
accuracy measurements, which shed some light on the influence of several tuning parameters and 
also show that the existing methods - although similar in spirit - exhibit remarkable differences 
in accuracy. We propose combinations of the individual components, mostly relying on the P 3 M 
approach, which we regard as most flexible. 



INTRODUCTION 

A challenging task in every computer simulation of par- 
ticles which are subject to periodic boundary conditions 
and long range interactions is the efficient calculation of 
quantities like the interparticle foregs or the interaction 
energies. The famous Ewald surroEl does a remarkable 
job in splitting the very slowly (not even uncondition- 
ally) converging sum over the Coulomb potential into 
two sums which converge exponentially fast. Still, this 
method suffers from two deficits: First, it is computation- 
ally demanding, since one part of the problem is solved 
in reciprocal space, thereby implying the need for several 
Fourier transformations. Second, the algorithm scales 
like A^ 2 (with N being the number of charged particles 
in the simulation box) or at best like N z / 2 , if one uses 
cutoffs whiph are optimized with respect to the splitting 
parametem. 

Several methods have been proposed to tackle the first 
problem, e.g. tabulation of the complete Ewald potential! 
or the use of polynomial approximations, in particular ex- 
pansion of the non-spherical contributions to the Ewald 
potential in cubic harmonicscl'B. Apart from the diffi- 
culty of a computational overhead which might strongly 
increase with the desired accuracy, all these methods do 
not solve the second problem: the unfavorable scaling 
with particle number. 

The essential idea is not to avoid the Fourier trans- 
forms but to modify the problem in a way that permits an 
employment of the Fast Fourier Transformation^ (FFT), 
thereby reducing the complexity of the reciprocal part of 
the Ewald sum to essentially order NlogN. If the real 
space cutoff is chosen small enough, this scaling applies to 
the complete Ewald sum. Since the FFT is a grid trans- 
formation, there are discretization problems to be solved 



and corresponding discretization errors to be minimized. 

At present there exist several mesh implementations 
of the Ewald sum - similar in spirit but different in de- 
tail. In this article we will focus on the original partic- 
le-particlejparticle-mesh (P 3 M) method of Hockney and 
EastwoodQ and two variants, namely, the particle mesh 
Ewald (PME) method of Dardqa et. ala and an upgrade 
of the latter by Essmann et. aln, which we will refer to 
as SPME (the "S" stands for "smooth"). 

There have been some uncertainties in the literature 
concerning the relative performance of these methods and 
it has been shown previouslytJ that the P 3 M approach 
- the oldest of the three - is actually the most accurate 
one and should be the preferred choice. However, since 
in this reference the PME method was combined with a 
disadvantageous charge assignment scheme and the more 
recent SPME could not be considered, we found it worth- 
while to test again these three methods under similar well 
posed and reproducible conditions and a larger number 
of tuning parameters. 

The original literature on particle mesh routines is 
mostly not easy to digest for the layman, obscured by 
the fact that the various authors approach the problem 
from different directions and use different notations. In 
this article we try to present a unified view of the common 
methods and analyze in detail the ingredients compris- 
ing them. By this we want to uncover the large number 
of possibilities for combining the different parts, thus al- 
lowing a judicious balance of accuracy, speed and ease 
of implementation. Moreover, we show that due to some 
subtle interdependencies not all combinations are advan- 
tageous, although they might appear promising at first 
sight. 

This paper is structured as follows: In the first section 
we briefly review the idea of the Ewald sum and pro- 
vide the most important formulas. Then we describe in 
some detail the steps which must be carried out if FFT 
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algorithms are to be employed for the Fourier transfor- 
mation, namely: charge assignment onto a mesh, solving 
Poisson's equation on that mesh, differentiating the po- 
tential to obtain the forces and interpolating the mesh 
based forces back to the particles. All these steps can 
be performed in different ways and in the following sec- 
tion we investigate in detail their accuracy, in particular, 
we compare P 3 M, PME and SPME with respect to their 
root mean square error in the force. Based on our the- 
oretical and numerical investigations, we find that the 
most accurate and versatile routine is the P 3 M method, 
supplemented by ingredients (dealing with the differenti- 
ation) of the other two approaches. 

The important task of an optimal tuning of the pa- 
rameters - especially the Ewald parameter a - is sat- 
isfactorily solved for the standard Ewald and the PME 
method, since there exist accurate analytic estimates for 
the root mean square error in the force. We will tackle 
the corresponding problem for P 3 M in a forthcoming 
publication!^. 



THE EWALD SUM 

There are many examples of long range interactions 
which can be treated by Ewald techniques, but in this 
article we will solely be concerned with Coulomb point 
charges, i.e. with an interaction potential 1/r. Consider 
therefore a system of N particles with charges qi at po- 
sitions Yi in an overall neutral and (for simplicity) cubic 
simulation box of length L and volume Vb = L 3 . If peri- 
odic boundary conditions are applied, the total electro- 
static energy of the box is given by 
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The sum over n takes into account the periodic images 
of the charges and the prime indicates that in the case 
i = j the term n = must be omitted. Of course r.y = 
i"i — r j , and our unit conventions are shortly described in 
Appendix [A]. 

Strictly speaking, since this sum is only conditionally 
convergent, its value is not well defined unless one spe- 
cifies the precise way in which the cluster of simulation 
boxes is supposed to -fill the R 3 , i.e. its shape (e.g. ap- 
proximately sphericaO) and the conditions outside the 
cluster (e.g. vacuum or some_dielectric). A thorough dis- 
cussion is given elsewhereo'Ej. 

The slowly decaying long range part of the Coulomb 
potential renders a straightforward summation of Eqn. 
([[]) impracticable. The trick is to split the problem into 
two parts by the following trivial identity: 



1 = f(T) , 1 - f(r) 
r r r 



(2) 



The underlying idea is to distribute the two main com- 
plications of the Coulomb potential - its rapid variation 



at small r and its slow decay at large r - between the 
two terms by a suitable choice of /. In particular: 

• The first part ^r- should be negligible (or even 
zero) beyond some cutoff r max , so that summation 
up to the cutoff is a good approximation to (or 
the exact result of) this contribution to the total 
electrostatic potential. 

• The second part 1 ~^ r ^ should be a slowly varying 
function for all r, so that its Fourier transform can 
be represented by only a few k- vectors with |k| < 
fcmax- This permits an efficient calculation of this 
contribution to the total electrostatic potential in 
reciprocal space. 

Since the field equations are linear, the sum of these 
two contributions gives the solution for the potential of 
the original problem. However, the two requirements on 
the function / mentioned above leave a large freedom of 
choicai-3. The traditional selection is the complementary 
error function erfc(r) :— 2tt~ 1 / 2 / r °° dt exp(— t 2 ), which 
results in the well known Ewald formula for the electro- 
static energy of the box: 



E = E {r) + E {k) + E (s) + E (d) 



(3) 



where the contribution from real space E^', the contri- 
bution from reciprocal space E^ k \ the self energy E^ 
and the dipole correction E^ are given by 
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and the Fourier transformed charge density p(k) is de- 
fined as 



p(k) = f d 3 r p(r) 
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The inverse length a, which we will refer to as the Ewald 
parameter, tunes the relative weight of the real space and 
the reciprocal space contribution, but the final result is of 
course independent of a. The k- vectors form the discrete 
set {2nn/L : n e Z 3 }. 

The form (0) given for the dipole correction assumes 
that the set of periodic replications of the simulation box 
tends in a spherical way towards an infinite cluster and 
that the rHiadium outside this sphere is a homogeneous 
diclcctricHil3 with dielectric constant e'. Note that the 
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case of a surrounding vacuum corresponds to e' = 1 and 
that the dipole correction vanishes for metallic boundary 
conditions, since then e' = oo. Note also that this term is 
independent of a, which again shows that it is not specific 
to the Ewald sum but more generally reflects the prob- 
lems inherent to the conditional convergence of the sum 
in Eqn. ([j]). Some complications regarding the cornet 
implementation of this term are discussed by CaillolEj. 

The advantage of rewriting Eqn. ([!]) this way is that 
the exponentially converging sums over m and k in (Q,^) 
allow the introduction of relatively small cutoffs with- 
out much loss in accuracy. Typically one chooses a large 
enough as to employ the minimum image convention in 
Eqn. (Q). It is important to realize that at given real- 
and reciprocal space cutoffs there exists an optimal a 
such that the accuracy of the approximated Ewald sum 
is as high as possible. This optimal value can be deter- 
mined easily with the help of the excellent estimates for 
the cutoff errors derived by Kolafa and PerranxJ - essen- 
tially by demanding that the real- and reciprocal space 
contribution to the error should be equal. 

Some more insight into the Ewald sum can be gained 
by the following considerations. Let g(k) := An/k 2 be 
the Fourier transformed Green function of the Coulomb 
potential 1/r and 7 (k) := exp(-fc 2 /4a 2 ). Then Eqn. (||) 
can be rewritten as follows: 



£fl(k)7(k)/5(k) t 
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Here <p ' (vj ) is the electrostatic potential at the point r j 
due to the second term in Eqn. (0) and by definition it is 
clear that its Fourier transform is given by 



(k) = ff(k) 7 (k)p(k) 



(10) 



As is known, products in reciprocal space correspond to 
convolutions in real space. Hence Eqn. (|l0|) shows that 
the reciprocal space contribution to the electrostatic po- 
tential is created by a charge distribution which is ob- 
tained from the original point charge distribution by a 
convolution with a "smearing function" 7(r). 

For the standard Ewald sum 7 (r) is a Gaussian, i.e. 
7(r) = a 3 7r~ 3 / 2 exp(— a 2 r 2 ), but this is merely a conse- 
quence of choosing the splitting function / in Eqn. (0) to 
be the complementary error function. In fact, an alterna- 
tive methodic to motivate the splitting which was done 
in Eqn. (0) is to replace the point charge distribution p by 
a screened charge distribution p — p* 7 and compensate 
this screening by adding the smeared charge distribution 
p*7- (The star denotes the convolution operation.) From 
a mathematical point of view these two interpretations 
are perfectly equivalent: Instead of splitting the potential 
one splits the charge density. 



At this point a word of caution seems appropriate: 
Whereas the electrostatic potential depends linearly on 
the charge density, the electrostatic energy does not. 
Thus, calculating the energies resulting from the charge 
densities p — p-kj and p* 7 and adding these contributions 
together would not give the energy of the charge density 
p. Consequently, E^ is not the electrostatic energy of a 
charge density p* 7 but the Fourier space contribution to 
the electrostatic energy of the charge density p. We want 
to make this subtle point more clear by writing down 
the energy explicitly. If we denote with (f> p the potential 
originating from p, we have due to the linear dependence 



of (j) p on p an equation like <p p = 



j p-p*i 
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Hence 



we can obtain for the electrostatic energy the following 
expression: 

E' = \f d 3 r p(r)<f, p (r) 



2 / dr P( Y ) 



/ dVp(r)^(r) + - / d J rp(r)^ 7 (r) (11) 

The two terms in the last line are the real space and the 
Fourier space contribution to the energy, but neither of 
them can be interpreted as the energy of a charge dis- 
tribution p — p * 7 or p * 7! Moreover, the quantity E' 
contains unphysical self energy contributions, i.e. energy 
due to the interaction of a charge (or 7-smeared charge) 
with itself. In the actual Ewald sum the self energy con- 
tribution of the real space part is canceled by omitting 
the term m = for i = j in Eqn. (|4|), whereas the self en- 
ergy contribution of E^ k > must be subtracted separately 
(this is the origin of the term E^). 

Finally, the force Fj on particle i is obtained by dif- 
ferentiating the electrostatic potential energy E with re- 
spect to rj, i.e. 
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Using Eqns. (|-|) one obtains the following Ewald for- 
mula for the forces: 



(13) 



with the real space, Fourier space and dipole contribu- 
tions respectively given by: 
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Since the self energy from Eqn. (g) is independent of 
particle positions, it does not contribute to the force. 

EWALD SUMMATION ON A GRID 

Performing the Fourier transformations inherent to the 
reciprocal space part of the Ewald sum by FFT routines 
is by no means a straightforward business. First, the 
point charges with their continuous coordinates have to 
be replaced by a grid based charge density, because the 
FFT is a discrete and finite Fourier transformation. Sec- 
ond, it is neither obvious nor true that the best grid ap- 
proximation to the continuum solution of Poisson's equa- 
tion is achieved by using the continuum Green function. 
Third, there are at least three ways for implementing 
the differentiation needed in Eqn. (p"2|), which differ in 
accuracy and speed. And fourth, the procedure of as- 
signing the forces calculated on the mesh back to the 
actual particles can - under certain circumstances - lead 
to unwanted violations of Newton's third law, which can 
be anything between harmless and disastrous. 

The four steps involved in a particle mesh calculation 
are sources for various kinds of errors, originating e.g. 
from discretization, interpolation or aliasing problems 
(with the latter we want to denote inaccuracies result- 
ing from the fact that a finite grid cannot represent arbi- 
trarily large k- vectors). Since these contributions are not 
independent of each other (reducing one might enhance 
another), the only reasonable demand is the minimiza- 
tion of the total error at given computational effort. 

One of our aims is to compile some of the possibilities 
for each step, in order to draw a comparison between the 
three mesh implementations mentioned in the introduc- 
tion - PME, SPME and P 3 M. Like the Ewald sum, all 
these algorithms can be extended to a triclinic simula- 
tion cell by reverting to general dual basis vectors and 
one can also use a different number of grid points along 
each direction. However, in order to keep the notation 
simple, we restrict to the case of a cubic box and employ 
the same number of mesh points in each direction. How 
the generalizationSj-jCan be done is described e.g. in the 
references on PMEEl or SPMEEJ. 



Charge Assignment 

The actual procedure of assigning the charges to the 
grid can be written down very easily. We will first discuss 
the one-dimensional case, i.e. particles with coordinates 
x € [0; L] C R have to be assigned to the mesh points 
x p € M = {p h : p = 0, . . . , Nm — 1}, where Nm is the 
number of mesh points and h := L/Nm is their spacing. 
To keep the notation simple, we will abstain from explicit- 
ly taking into account that any x-value, which is outside 
[0;L], has to be folded back into this interval in order 
to conform to periodic boundary conditions. Rather, we 



assume that this is done as necessary, i.e. all calculations 
are to be understood "modulo i" . 

Define the even function W(x) such that the fraction 
of charge which is assigned to the mesh point x p due to 
a unit charge at position x is given by W(x — x p ). If the 
charge density of the system is p{x), then the mesh based 
charge density pm, defined at the mesh points x p , can be 
written as the following convolution: 

Pm(x p ) = i J dx W(x p - x) p(x) (17) 

The prefactor 1/h merely ensures that pm is in fact a 
density. Henceforth we will refer to any such W as a 
charge assignment function. 

The important question is: What properties should 
W(x) have in order to be a suitable choice? The fol- 
lowing wish list summarizes some desirable features: 

• Charge conservation, i.e. the fractional charges of 
one particle, which have been distributed to the 
surrounding grid points, sum up to the total charge 
of that particle. 

• Finite and if possible small support, because the 
computational cost increases with the number of 
mesh points among which the charge of each par- 
ticle is distributed. (The support of a real-valued 
function / defined on X is (the closure of) the set 
{x € X : f{x) ^ 0}, i.e. basically the range of 
values for which the function is nonzero.) 

• Localization of discretization errors, i.e. inaccura- 
cies in the force between two particles due to the 
discretization should become small with increasing 
particle separation. 

• Large degree of smoothness, i.e. the fractional 
charge of particle i which is assigned to some mesh 
point x p should be a smoothly varying function of 
the position of particle i. 

• Minimization of aliasing errors, i.e. since on a finite 
grid there is only space for a limited number of 
k- vectors, the charge assignment function should 
decay sufficiently rapidly in Fourier space. 

• Easy and transparent implementation. 

It is important to realize that these characteristics can- 
not be achieved all at the same time. Although some 
properties are positively correlated (e.g., a large degree 
of smoothness implies a fast decay in reciprocal space and 
thus minimizes aliasing errors), some other properties ex- 
clude each other (e.g., minimization of aliasing errors im- 
plies a sufficient localization in reciprocal space which is 
incompatible with a small support in real space). Thus, a 
good charge assignment function is always a compromise 
between these different demands. 
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FIG. 1. Thi]°d-»rder charge assignment function few the La- 
grange schemeB'EZI (solid line) and the spline schemed (dotted 
line, see also Eqn. ([[9|)). Both assignment functions have 
support [— | ft; | ft] and are piecewise quadratic. While for all 
odd assignment orders the Lagrange assignment function is 
discontinuous (for even assignment order it is continuous but 
not differentiable), the spline assignment function is in general 
P — 2 times continuously differentiable by construction. Note, 
however, that the Lagrange assignment function is optimized 
with respect to a different property, see Appendix [ 



We also want to stress that the choice of W(x) is not 
independent of the other decisions made for the mesh im- 
plementation. In Appendix |b] we show that if one sticks 
to the continuum Coulomb Green function in the mesh 
calculation, the requirement of localization of discretiza- 
tion errors is enough to restrict the charge assignment 
to a Lagrange interpolation yS^heme. This is in fact the 
combination used for PMEBtLZl (see also Appendix |b|). 
Hence, other influence functions can only be competitive 
if the Coulomb Green function is somehow adjusted at 
the same time. 

The choice for the charge assignment function of Hock- 
ney and Eastwood is as follows: In a P th order assign- 
ment scheme (i.e. the charge of one particle is distributed 
between its P nearest mesh points) define the Fourier 
transformed charge assignment function as 



kh/2 J 

Transforming this back to real space gives 



(18) 



(19) 



f-fold convolution 



with X[-| i] being the characteristic function of the in- 
terval [—5, 5], i.e. the function that is 1 within this in- 
terval and outside. Thus, e.g. by the central limit 
theorem, the P th order charge assignment function re- 
sembles (in this case) for increasing P more and more 



closely a centered Gaussian (with a variance P times as 
large as the variance of W^), but it has finite support 
[— Tfj =jf]. This assignment function is very smooth for 
large P, since it is .-a. spline of order P and thus P — 2 
times differentiablella. As a matter of convenience we 
decided to tabulate the corresponding charge fractions 
W p (p) ( x ) = W (P) (x- x p ) for P e {1, . . . , 7} in Appendix 

^ • n 

In Fig. HI the third order assignment functions for the 
Lagrange and the spline interpolation scheme are plot- 
ted. Note that while the spline function is in general 
P — 2 times continuously differentiable, the Lagrange as- 
signment function is not as smooth: Generally, for even 
assignment orders it is continuous, but the derivative is 
not, while for odd assignment orders it is discontinuous 
right away. Incidentally, for P = 1 and P — 2 both 
schemes coincide. 

The SPME method uses in essence the same charge 
assignment functions as the P 3 M-method, but this is dis- 
cussed more appropriately in the next section. 

Charge assignment in more than one dimension can be 
achieved by a simple factorization approach. E.g., the 
three-dimensional charge assignment function W(r) can 
be written as 

W(r) = W(x)W(y)W(z) (20) 

This is certainly not the only possibility^, but it is com- 
putationally advantageous. 

The generalization of Eqn. ( |l7| ) to three dimensions 
can be written as 
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(21) 
(22) 



In the last equation the reader should not confuse the 
coordinate of particle i, r,;, with the coordinate of mesh 
point p, r p . 



Solving Poisson's equation 

For the standard Ewald sum the Fourier space contri- 
bution to the electrostatic energy is given by Eqn. (|5|). 
How is this equation to be modified, now that we are 
working on a discrete mesh? r-. 

The simplest approach is used in the PME methods, 
where it is assumed that this equation is appropriate in 
the discrete case as well. The only difference is that the 
Fourier transformed charge density p from Eqn. (B) is re- 
placed by the finite Fourier transform of the mesh based 
charge density, /3m, which we define as 



/3m (k) := ft 3 ^2 Pm(t p ) e 



— i k-i 



(23) 
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where J^r em ^ s the sum over the (three-dimensional) 
mesh in real space and the k- vectors are from the corre- 
sponding Fourier space mesh. Of course, the way back 
to real space is also done by a (inverse) finite Fourier 
transform. (In order to distinguish between the usual 
and the finite Fourier transform, we indicate the latter 
by a hat and not by a tilde.) As discussed in the previous 
subsection and Appendix pi the usage of the continuum 
Coulomb Green function is best accompanied by a La- 
grange interpolation scheme for the charge assignment. 
See e.g. Petersentll for a tabulation of the corresponding 
polynomials and their implementation. 

A second algorithm, the SPME method, was presented 
by Essmann et. alxi. It uses a smooth charge assignment 
scheme and hence an adjusted Green function. The rea- 
soning is as follows: Starting with Eqns. (|],||) it is ar- 
gued that charge assignment onto the mesh can equiva- 
lently be viewed as interpolating exponentials of the form 
exp(ikx) at discrete grid points. This problem has a par- 
ticularly, elegant solution, the so-called exponential Euler 
splinescB. If x is the continuous particle coordinate, we 
have for even P (a recipe for the treatment of odd P can 
be found in the original SPME refercnceo) : 
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The function Mi p ^ is a cardinal-B-spline of order P, and 
Essmann et. alu give a recursive definition. They also 
point out that (in our notation for h — 1) the function 
Af( p ) is identical to the probability distribution of the 
sum of P independent random variables, each distributed 
uniformly on the unit interval. Since this distribution 
is given by the P-fold convolution of the characteristic 
function X[0;i] with itself, we can see by comparison with 
Eqn. ( |l9| ) that the charge assignment functions from the 
SPME and the P 3 M method are in fact identical up to 
a translation: M^ p \x) = W^ p \x - ^). Of course, 
A/( p ) has finite support [0; Ph], so the sum in Eqn. (p3) 
is actually finite. 

Note that is not an even function, and in the 

original reference on SPMEo the charge assignment dif- 
fers slightly from our Eqn. fl2^). However, the only effect 
of the translation is that the original system is repre- 
sented by a shifted mesh system, which is from a practi- 
cal point of view irrelevant, because this shift is undone 
in the back-interpolation (if accomplished with the same 
assignment function) . 

Now the following approximation for E^ can be de- 
rived by inserting Eqns. ( p4| , p5| ) into Eqn. (|J): 



h 3 p M (r p )[p M *G](r p ) 



(26) 



Here the star denotes the finite convolution 

[PM * G] (r p ) = h 3 J2 PM(r q )G(r p - r q ) (27) 
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(again, the periodic closure is not written down expli- 
citly) and the function G is given by its finite Fourier 
transform 



G(k) = B(k) 
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(28) 



with B(k) := \b(k x )b(k y )b(k z )\ 2 . Following Hockney and 
Eastwood we will refer to G as the influence function. 
The nice thing about G is that it is by construction in- 
dependent of particle coordinates and can therefore be 
precomputed. 

Eqn. ( |26] ) can be made plausible in the following way: 
G plays the role of a Coulomb Green function which 
has incorporated the "smearing" with the Gaussian 7. 
Hence, its convolution with the mesh based point charge 
density gives the mesh based electrostatic potential of 7- 
smeared charges. Multiplying this with the mesh based 
charge h 3 pM and summing over all mesh points gives the 
Fourier space contribution to the electrostatic energy up 
to a factor 1/2, which merely cancels some double count- 
ing. This should be compared to Eqn. (^) or the second 
term in Eqn. (pi]). 

As pointed out before, a charge assignment different 
from the Lagrange interpolation scheme can only be com- 
petitive, if the Coulomb Green function is changed at the 
same time. The replacement of the usual (and smeared) 
Green function 0*7 with the influence function G, which 
essentially differs by the additional prefactor B in Fourier 
space, achieves exactly that. Conversely, PME uses a La- 
grange interpolation scheme together with the unchanged 
Coulomb Green function, i.e. Eqns. ( f26| , p8|) with B = 1. 

Finally, the alias sum occurring in Eon. ( |28| ) is sub- 
stituted according to the following rulea: If Nm is the 
number of mesh points in each direction and the vector 
k on the left hand side is given by k = 27rn/i, n £ 
{0, ...,N M - I} 3 , then define G(k) = B(k)g(k')^(k'), 
where k' = 27rn'/L and n\ = rii for < i%i < Nyi/2 and 
n[ = ni — Nm otherwise (i — x, y, z). 

A third possibility - the so-called P-?M method - was 
presented by Hockney and EastwoodP: Their objective 
was an optimization of the influence function G in Eqn. 
(p6|), which causes the final result of the mesh calcula- 
tion to be as close as possible to the original continuum 
problem. So in order to proceed one first has to make 
the statement "as close as possible" more quantitative, 
and this can be done as follows: 

Take two particles with coordinates 17 and Y2 and de- 
fine r := 17 — Y2- The true force between these particles 
should be a function of r only, but in any mesh imple- 
mentation the actual force also depends on the positions 
of the particles relative to the mesh, say, on the position 



G 



of the first particle within its mesh cell (i.e. the origi- 
nal translational symmetry is broken by the mesh) . This 
suggests the following measure for the error: integrate 
the square of the difference between the calculated force 
F and the true reference force R over all values of r and 
average this quantity over all positions of e.g. the first 
particle within one particular mesh cell: 



Q~l f d 3 ri f d 3 r[F(r;n)-R(r)] s 

*C JVa JV b 



(29) 



Here V c = h 3 is the volume of one mesh cell. The solution 
of Poisson's equation is accomplished in essence by Eqn. 
( p6| ) (it is only written down somewhat differently), and 
the derivative in Eqn. ( |l2| ) is performed by applying finite 
difference operators to the mesh based electrostatic po- 
tential (see below). Since the discretization error Q can 
be regarded as a functional of G, the optimal influence 
function G pt can be obtained by setting the functional 
derivative of Q with respect to G to zero, i.e. 



SQ 
SG 



= 
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G — Gopt 



Starting from this idea, Hockney and Eastwood were able 
to derive the following expression for G op tB: 



Gopt(k) 



D(k) • E meZ 3 U 2 (k + f m)R(k + f m 



D(k)| 2 E meZ 3[/ 2 (k+» 



(31) 



Here D(k) is the Fourier transform of the employed dif- 
ferentiation operator (see next section and Appendix ^) , 
U(k) — W(k)/V c is the Fourier transform of the charge 
assignment function divided by the volume of one mesh 
cell and R(k) is the Fourier transform of the true refer- 
ence force, given by 



R(k) = -<kfl(k)7(k) 



(32) 



Note that this differs from ihe expression from the book 
of Hockney and EastwooaO, who use 7 2 instead of 7. 
The reason is that Eqn. (|3^) describes the true reference 
force between a 7-smeared charge and a point charge, 
while Hockney and Eastwood choose a slightly different 
approach in which they need the force between two 7- 
smeared charges. Also, we keep the factor 4tt in the 
Fourier transformed Green function g and do not to hide 
it somewhere else. 

The alias sums over m in Eqn. ([u|) are typically well 
converged for |m| < 2 and the sum in the denominator 
could even be done analytically. Again we want to em- 
phasize that the calculation of the influence function has 
to be done only once prior to the actual simulation and 
thus does not produce any runtime overhead. Note also 
that the expression ( p8| ) differs from the optimal form 
( [31] ) and hence cannot be optimal. 



A final word concerning the implementation: Although 
the convolution pm * G in Eqn. ( pfj| ) is a nice and com- 
pact notation, the whole purpose of these particle mesh 
routines is to employ the convolution theorem and use 
efficient FFT routines to calculate pm * G. The central 
steps are thus: 

• Calculate the finite Fourier transform /3m of the 
mesh based charge density pm- 

• Multiply pm with the precomputed Fourier space 
representation of the influence function, G. 

• Apply an inverse finite Fourier transform to this 
product to end up with the finite convolution of 
Pm with G. Formally this can be symbolized as 



p M * G = FFT 



FFT [p M ] x FFT [G] 



(33) 



Note that in this way one only needs G to calculate 
Pm * G but actually never G itself. 

This is the important part which all particle mesh al- 
gorithms have in common. The various methods differ 
e.g. in their choice of G, the assignment function W or 
the implementation of the derivative in Eqn. (p^). 



Differentiation 

After the calculation of the electrostatic energy, the 
forces on the particles are obtained by differentiation ac- 
cording to Eqn. (|l2]). However, for the Fourier space part 
of particle mesh methods there are several possibilities 
to implement this procedure. In other words: there exist 
several possible substitutes for Eqn. (|l5|), in particular 

1. Differentiation in Fourier space. 

2. Analytic differentiation of the assignment function 
in real space. 

3. Discrete differentiation on the mesh in real space. 

Differentiation in Fourier space is easy, since it merely 
involves a multiplication with the Fourier transformed 
differentiation operator D(k), which is a fast, local and 
accurate operation. Although one might want to use 
Fourier transforms of discrete difference operators - al- 
lowing for the fact that one is actually working on a 
mesh - the best results are obtained when the Fourier 
transform of the usual differential operator, namely ik, 
is employed. Therefore we will refer to this method as 
ik-differentiation. The basic idea is not to calculate the 
mesh based electrostatic potential (j^ k \v p ) via Eqn. ( |33| ) 
but the mesh based electric field E(r p ) by the following 
simple change to this equation: 



7 



dr p 
FFT 



dr p 

ik x p M x G 



(34) 



This method is employed in the PME algorithm and, as 
shown later, leads to the most accurate force calculations, 
if it is used in conjunction with the optimal influence 
function from Eqn. (|3l"|). Note, however, that since k is 
a vector, there are in fact three inverse three-dimensional 
Fourier transforms to be calculated in Eqn. (|34|), which 
is obviously computationally demanding. 

The electrostatic energy calculated on the mesh de- 
pends on the particle coordinates through the arguments 
of the charge assignment function W. As the creators of 
the SPME method point outB, a smooth charge assign- 
ment scheme permits an analytic differentiation of the 
energy, since the quantity pm, which contains the parti- 
cle coordinates i v depe nds in a differentiable way on the 
Yi. Using Eqns. (|l^ , p6| ) and the fact that G is indepen- 
dent of particle coordinates and an even function (since 
G is even), one can derive 



d_l 



^2 h 3 p M {r p )[pM*G](r p 



r p ef 



h <L [-^r( r p)PM(r q )+pM(r P )-^-(r q )) x 



2'" ^ V dr. 
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X 



xG(r p - r 9 ) 

I 

dr 

\{g{v p - r q ) + G(r q - r p ) 

P_ 
dr 



- h3 H ^p)[PM*G](r p ) 



(35) 



From Eqn. ( pi|) it is obvious that the array [dpM/ dri]{r p ) 
is essentially obtained by a charge assignment scheme 
which uses the gradient of the assignment function W and 
can thus be calculated conveniently at the same time as 
Pm- Since only one Fourier transform back to real space 
is necessary, this procedure is indeed very fast. Unfortu- 
nately, this differentiation scheme leads to a small ran- 
dom particle drift, since momentum is not conserved any 
more (see next section). Although the total momentum 
of the simulation box can be kept constant by subtract- 
ing the mean force jj from each particle, the small 
reduction in the accuracy of the particle forces due to 
these local random fluctuations can only be compensated 
marginally by this global correction. 

A third possibility for implementing the derivative in 
Eqn. jl^ ) is the use of finite difference operators, which 
calculate the force on one mesh point from the poten- 
tial at the neighboring mesh points. This is basically 
the method which is favored by Hockney and Eastwood 
for P 3 M. Higher accuracy is achieved by considering not 



only the nearest neighbors but also mesh points farther 
away, i.e. using linear combinations of nearest neighbor, 
next nearest neighbor etc. difference operators. In Ap- 
pendix |c| we show, how these approximations are con- 
structed systematically. In the P 3 M-method the Fourier 
transforms of these operators are needed for the calcu- 
lation of the optimal influence function (|3l]). This ap- 
proach as well needs only one Fourier transformation 
back to real space, like in the method of analytic differen- 
tiation. But unlike the latter, it conserves momentum (if 
the difference operators are chosen correctlyO) and thus 
has no problems with spurious particle drifts and result- 
ing errors in the force. However, using the neighboring 
points is a nonlocal approach and increasing its accuracy 
can only be done by taking into account more neighbors 
- which makes it even more nonlocal and more costly. 

Obviously, there is no unique optimal way for doing 
the differentiation. Each approach joins together advan- 
tages and drawbacks which have to be balanced against 
each other under the constraint of required accuracy and 
available computational resources. Let us make just one 
example: If the required accuracy is not very high, using 
only nearest neighbors for the discrete differentiation on 
the mesh might be accurate enough. Certainly, multipli- 
cation in Fourier space by ik gives better results, but let 
us assume that this approach is actually slower due to the 
two additional Fourier transformations. However, if the 
required accuracy increases, the finite difference approx- 
imation calls for more neighbors and thus becomes more 
and more costly, whereas the ik-approach right away 
gives the best result possible by discrete differentiation. 
This is because increasing the order of the differentiation 
scheme means that in Fourier space the transformed op- 
erators approximate ik to higher and higher truncation 
order (actually, that is how these approximations are con- 
structed, see Appendix |cj) . In other words, accepting the 
two additional Fourier transformations can be competi- 
tive. Moreover: The method of analytic differentiation 
could be faster than the discrete difference method even 
for J = 1. Thus, in cases where the latter is less accurate 
than analytic differentiation, there is no reason for using 
it. 

Whether there exists a break-even point between these 
methods and - if yes - where it is located can depend on 
the tuning parameters like mesh size and interpolation 
order as well as on the details of the implementation or 
the computational facilities one is working with. A gen- 
eral statement seems to be difficult. 



Back-interpolation 

At some stage of any particle mesh method a back- 
interpolation of the mesh based results to the actual par- 
ticles is necessary. As we have seen in the last subsection, 
this can be done before or after the Fourier transforma- 
tion back to real space, i.e. the mesh points can con- 
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tain either the potential or the components of the electric 
field. 

Basically, this back-interpolation is done in a similar 
way as the distribution of the charges to the mesh at 
the beginning of the calculation: via some assignment 
function W. E.g., the force on particle i is given by 

F, =QiJ2 E ( r p)^(r* - t p ) (36) 

r p GM 

with E(r B ) being the electric field on mesh point r p from 
Eqn. (p^). The interpretation of Eqn. (|3^ ) is the follow- 
ing: Due to the discretization each particle is replaced 
by several "sub-particles" , which are located at the sur- 
rounding mesh points and carry a certain fraction of the 
charge of the original particle. The force on each sub- 
particle is given by its charge times the electric field at 
its mesh point, and the force on the original particle is 
the sum of the forces of its sub-particles. 

From a technical point of view it is convenient to use 
the same function W for the assignment onto and from 
the mesh, because if in the first step one does not only cal- 
culate the total charge accumulated at some mesh point 
but additionally memorizes, to what extend the individ- 
ual particles contributed to this charge, the interpolation 
back can be done without a single function call to W. 

However, there is also a more subtle reason which sug- 
gests a symmetric interpolation, and this is related to the 
conservation of momentum. As demonstrated by Hock- 
ney and EastwoodQ, the force which a particle acts onto 
itself is zero and Newton's third law is obeyed (up to 
machine precision), if 

• charge assignment and force interpolation are done 
by the same function W and 

• the approximations to the derivatives are correctly 
space centered. 

The second requirement states that if the electric field 
at some mesh point r p can formally be written as 
Er„ d ( r p> r q)PM(r q ), then d(r p , r q ) = -d(r g , v p ). 

The method of analytical differentiation mentioned in 
the last subsection does not use mesh based derivatives 
and this is the approach chosen in the implementation 
of SPMEcl. In their paper the authors state that the 
sum of the electrostatic forces on the atoms is not zero 
but a random quantity of the order of the rms error in 
the force. We believe that these fluctuations have their 
origin in a violation of the above conditions, although 
strictly speaking these are only sufficient conditions for 
momentum conservation. 

For a more detailed discussion of related effects and the 
connection between momentum conserving and energy 
conserving methods see Hockney and EastwoodQ. 



INVESTIGATING THE ACCURACY 

An investigation of the errors connected with particle 
mesh Ewald methods is important for several reasons. 
First, the complete procedure of discretization introduces 
new sources of errors in addition to the ones originating 
from real and reciprocal space cutoffs. Second, compar- 
ing the efficiency of different mesh methods is only fair if 
it is done at the same level of accuracy. And third, the 
tuning parameters should be chosen in such a way as to 
run the algorithm at its optimal operation point. 

However, there is no unique or optimal measure of ac- 
curacy. If molecular dynamics simulation are performed, 
the main interest lies in errors connected with the force, 
while in Monte Carlo simulations one is concerned with 
accurate energies. In the simulation of ensemble aver- 
ages it is the global accuracy - measured e.g. by root 
mean square quantities - which is important, but in the 
simulation of rare events local accuracy and maximal er- 
rors are also relevant. Errors in the force can be due to 
their magnitude or due to their direction. And finally, 
one might be interested in absolute or relative errors. 

Whatever quantity one decides to look at, it can be 
investigated as a function of system parameters like par- 
ticle separation or distribution, tuning parameters like 
a, mesh size or interpolation order and components of 
the algorithm, e.g. interpolation or differentiation scheme 
or splitting function f(r). Obviously this gives rise to 
a very large number of combinations. In other words: 
The corresponding parameter space is large and nontriv- 
ial, i.e. general statements concerning the performance 
of one method can usually not be extracted from low- 
dimensional cuts through this space, because different 
methods scale differently with respect to their parame- 
ters. 

Nevertheless, we want to present some numerical accu- 
racy measurements at important points of this parameter 
space for the following reasons: As we pointed out, there 
are several options for the implementation of each step of 
a mesh calculation - e.g. three ways for doing the deriva- 
tive in Eqn. (p^). This freedom of choice and its impact 
on the overall accuracy has not been systematically in- 
vestigated so far, although a qualitative understanding 
of at least typical influences of the different parts on the 
performance permits a judicious assessment and compar- 
ison of the resulting algorithms, in particular P 3 M, PME 
and SPME. We want to show which combinations are 
attractive and which should definitely be avoided. And 
finally we want to present easily reproducible measure- 
ments which should allow the reader a comparison with 
his own implementations of particle mesh Ewald routines. 
However, we will not present large accurate tables, which 
provide an easy way for tuning these algorithms under 
all circumstances. On the contrary, we want to encour- 
age any potential user to perform some of these simple 
measurements on his own and thereby not only gaining 
insight but also the possibility to optimize his tuning pa- 
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rameters. We want to stress that parameters which are 
only roughly estimated or even historically handed down 
should be taken with great care. 



One possible measure of accuracy 

In this article we will solely be concerned with one 
measure of accuracy, namely the root mean square (rms) 
error in the force, given by 



AF := 



\ 



1 



N 

1=1 



(37) 



Where F$ is the force on particle number i calculated 
via some mesh method and F| xa is the exact force on 
that particle, calculable e.g. by a well converged stan- 
dard Ewald sum. There exist error estimates for the real 
space and Fourier space contribution to this errox-ior the 
standard Ewald sumli3 and for the PME methocO which 
greatly simplify the determination of the optimal value 
of a. 



Error as a function of a 

We investigated the rms error ( |37|) for a system of 
100 particles (50 carry a positive and 50 a negative unit 
charge) , which were randomly placed within a simulation 
box of length L = 10, as a function of the Ewald param- 
eter a. In order to make our results fully reproducible, 
we describe in Appendix [D], how our actual random con- 
figuration was generated. 

For small a the result of the Ewald sum (or any of 
the described particle mesh methods) is dominated by 
the real space contribution (ji|) while for large a it is the 
Fourier part (|^) which is important. This is a simple 
consequence of the fact that in the real space sum a oc- 
curs in the numerator of the exponential function (or - to 
be precise - of the complementary error function) while 
in the Fourier space sum it occurs in the denominator 
and thus influences the decay of both contributions in 
a converse way. Hence, at given cutoffs, the same ap- 
plies to the errors. Since with increasing a the real space 
contribution becomes more accurate while the Fourier 
space contribution degrades in accuracy, one can expect 
an optimal a to exist at which the total error is minimal. 
This is approximately at the point where real space and 
Fourier space errors are equal. Since the different mesh 
methods we investigate all coincide in the treatment of 
the real space part, their errors should all be the same 
for sufficiently small a. 

In Fig. H we plot the rms error of the force as a function 
of a, which was obtained by investigating our system 
with various mesh methods. They all share a mesh size 
of Am = 32 (and thus have 32 3 mesh points in total), 
an interpolation order P = 7 and a real space cutoff 



r max = 4. We find indeed the general features described 
above, like a low accuracy for very small or very large 
values of a and an optimal value in-between. However, 
the various methods differ considerably in their accuracy. 
(Note that in this and the following figures the vertical 
scale is logarithmic!) 

The solid line corresponds to PME. This method com- 
prises some elements which make one think about pos- 
sible improvements: a not very smooth charge assign- 
ment scheme (namely, the Lagrange interpolation) and 
the use of the plain continuum Green function. There is 
clearly no obvious advantageous replacement for the lat- 
ter, but it is easy to replace the Lagrange scheme by the 
smooth spline interpolation (by just changing the assign- 
ment function). Yet, the result of this supposed improve- 
ment, shown in line 2, is in fact disappointing. As we 
mentioned several times, the continuum Green function 
is best accompanied by a Lagrange interpolation scheme, 
because this leads to a cancellation of certain discretiza- 
tion errors. Changing the assignment scheme destroys 
this effect and the resulting error shatters the desired 
improvement in accuracy completely. 

Upgrading PME requires a proper treatment of both 
elements - charge assignment and Green function. This 
is in fact what the remaining two algorithms (SPME and 
P 3 M) accomplish. Since they both use a smooth spline 
interpolation, they are both potential candidates for an- 
alytic differentiation. In fact, the SPME method, as de- 
scribed in the original publication^, chooses this imple- 
mentation of the derivative, because it is very fast (line 
3). Nevertheless, the ik- method is still possible and leads 
to an even better result (line 4) , which admittedly has to 
be payed with two additional FFT calls. Analytically 
differentiated P 3 M gives an error almost identical to an- 
alytically differentiated SPME, but if one implements the 
ik-derivative, P 3 M improves a little bit on SPME. From 
a theoretical point of view the latter is not too much 
surprising: After all, if P 3 M uses an optimal differen- 
tiation (in view of accuracy) and an optimal influence 
function, it can be expected to constitute a kind of lower 
bound for the error. However, if the optimal differenti- 
ation is replaced by the analytic differentiation, a new 
source of error appears (namely, the random force fluc- 
tuations described in the section on back-interpolation). 
If this contribution dominates, the fact that P 3 M uses 
a better influence function than SPME cannot make a 
large difference. In our case the analytically differen- 
tiated SPME is a factor 9.2 more accurate than PME, 
while the ik-diffcrcntiated P 3 M method is more accu- 
rate than PME by a factor of about 33. However, one 
must realize that SPME and P 3 M have different execu- 
tion times, since P 3 M needs two additional FFT calls 
compared to SPME. But apart from the analytically dif- 
ferentiated curves all methods summarized in Fig. |2j need 
exactly the same time for a mesh calculation. This comes 
from the fact that the methods differ only in parts which 
normally are tabulated anyway, like the influence func- 
tion. 
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FIG. 2. Comparison of different mesh methods: The rms error AF from Eqn. @ for a system of 100 charged particles 
randomly distributed within a cubic box of length L — 10 (see Appendix ^) is shown as a function of the Ewald parameter a 
for 6 mesh algorithms, which all share Nm = 32, P — 7 and r m ax = 4. Line 1 is PME. Line 2 corresponds to an algorithm 
which is obtained from PME by retaining the continuum Green function but changing to the spline charge assignment. Lines 
3 and 4 are analytically and ik-differentiated SPME respectively and line 5 and 6 are analytically and ik-differentiated P 3 M 
respectively. Note the logarithmic vertical scale in this and the following figures. 



There is another surprising thing to note about SPME: 
For the chosen values of Nm and P the curves for PME 
and (analytically differentiated) SPME intersect, i.e. the 
latter is not necessarily more accurate. It could be argued 
that at least for the optimal value of a SPME is better, 
but this optimal a of course depends on the real space 
cutoff r max as well. If this cutoff is decreased, the real 
space contribution to the error is increased. In fact, using 
the estimate of Kolafa and Perramtj one finds that at the 
intersection point of PME and SPME this contribution 
will have the same size as the Fourier space contribution 
for (in this case) r max 1. Thus, for even smaller val- 
ues of r max PME would actually be more accurate than 
SPME. 

Now that we have compared various particle mesh 
methods, we want to examine in a little more detail some 
parts of the algorithm. We will always use the P 3 M 
method for illustration. Corresponding plots for PME or 
SPME would look qualitatively very similar and hence 
are not presented in the sequel. 

In Figure || we took all parameters of the ik- 



differentiated P 3 M method from Fig. || but varied the 
charge assignment order from P = 1 to P = 7. Increas- 
ing P improves the accuracy by more than three orders of 
magnitude (from P = 1 to P = 7). However, the reward 
for going from P to P + 1 is larger for small P. Note also 
that the optimal value of a depends on P. 

In Figure ^ we fix the order of the charge assignment 
scheme to P = 3 and vary the number Nm of Fourier 
mesh points. Qualitatively the behavior is similar to Fig. 
^ Improving the method reduces the error and shifts 
the optimal a to the right. Note that from a computa- 
tional point of view Figs. || and |] are sort of conjugate: 
The accuracy depends on both Nm and P, but increas- 
ing one parameter does not influence the performance 
of the other. In other words: The charge assignment 
scales as P 3 independent of Nm and the FFT scales as 
(./Vm log-ZViyi) 3 independent of P. Optimal performance 
requires a suitable combination of Nm and P. 

Next we investigated the differentiation scheme. To 
this end we employed the P 3 M method with Nm = 32 
and P — 7 and used various orders J of the mesh based 
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FIG. 3. Influence of the charge assignment order: The rms 
error AF for our model system from Appendix^ is calculated 
for the ik-differentiated P 3 M method with = 32 and 

r max = 4. From top to bottom the order P of the (spline) 
charge assignment scheme is increased from 1 to 7. 



FIG. 5. Influence of the differentiation scheme: The rms er- 
ror AF for our model system from Appendix [5] is calculated 
for the P 3 M method with N M = 32, P = 7 and r max = 4. 
Shown are 6 mesh based approximations to the differenti- 
ation operator ik (from top to bottom: A' 1 - 1 , . . . , A' 6 ' , see 
Appendix |c|) as well as the result for ik itself (lowest curve, 
solid line). 




1.5 

FIG. 4. Influence of the mesh size: The rms error AF 
for our model system from Appendix |5] is calculated for the 
ik-differentiated P 3 M method with P — 3 and r maK — 4. 
From top to bottom the mesh size is given by 4,8,16,32,64 
and 128. (Note that the total number of mesh points in this 
three-dimensional system is given by (A^m) 3 -) 



approximation to the difference operator (see Appendix 
^jj). (Actually, the calculations were done by a multi- 
plication in Fourier space with the transformed approx- 
imations D^' 7 ).) The result is shown in Fig. |[ which 
looks pretty much like Fig. [| but was generated quite 
differently. With increasing order of the difference ap- 
proximation the errors decrease. However, the result of 
the ik-diffcrentiation scheme forms a lower bound to the 
error of this method. (After all, ik is the Fourier rep- 



FIG. 6. Comparison of a mesh method with the standard 
Ewald sum: The rms error AF for the Ewald (line 1) and 
PME (line 2) method are calculated for our model system 
from Appendix |d| The parameters for PME are the same as 
in Fig. H and the Fourier space cutoff for the Ewald sum was 
set to femax = 20 x 2tt/L. This value is interesting to compare 
with the PME method, because it corresponds to the same 
number of k-vectors (since |7r2B 3 i « 32 3 ). Also shown is the 
estimate for the real space errortJ (line 3), the Fourier space 
error for EwakU (line 4, we used the slightly better estimate 
from PeterserJiZI) and the Fourier space error for PMEliJ (line 
5). Note that the estimates for the Ewald sum can hardly be 
distinguished from line 1. 



resentation of the exact differential operator, and in the 
standard Ewald sum the differentiation is also done this 
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way, compare Eqn. (p"5|).) Here the bound is reached in 
the minimum at J = 7, so further improving the dif- 
ferentiation order is of no use at all. Of course, if the 
accuracy of the lower bound is smaller (e.g., because the 
charge assignment order is lower) the zk-bound will be 
reached already by smaller values of J. Note that in this 
example the method of analytic differentiation gives ap- 
proximately the same accuracy as a fifth order difference 
scheme (compare to Fig. |2|). Since analytic differentiation 
is much faster, it should be preferred to the finite differ- 
ence approach in cases where the latter is less accurate 
anyway. 

The last part of this section deals with the determi- 
nation of the optimal a-value. There exist rather good 
estimates for the real-, and reciprocal space error of the 
standard Ewald sumll3 and the reciprocal space error of 
the PME methodlJ. The optimal a-value of these two 
methods and the corresponding accuracy can be obtained 
very precisely by just calculating the intersection point 
of the real- and corresponding reciprocal space estimates. 
Their high quality is clearly demonstrated in Fig. ^. The 
existence of these formulas is certainly a big advantage 
of the PME method, since it permits an a priori deter- 
mination of the optimal operation point as a function of 
system specifications (like box length, particle number or 
valence) or method parameters (like mesh size or assign- 
ment order) . The elaboration of a similar error estimate 
for the P 3 M method is currently pursued and will be pre- 
sented in a forthcoming publicationt 1 ! This is basically 
the last step which is missing to advocate P 3 M as the 
most accurate and versatile Ewald mesh method. 

A final word concerning the accuracy of mesh meth- 
ods compared to the Ewald sum: The optimal a value 
for a standard Ewald summation of our system with 
k m ax = 20 x 2ir / L (which thus has the same number of 
k-vectors, because |tt20 3 » 32 3 ) is approximately 1.25 
and the corresponding total error is of the order 5 x 10~ 12 
(see Fig. |J). Although much optimization effort has been 
put into mesh methods in order to reduce errors, we must 
face the fact that one generally loses many orders of mag- 
nitude in accuracy due to discretization. So if high ac- 
curacy is essential but speed is not an issue, the conven- 
tional Ewald method is unsurpassed: it is much easier to 
program and the desired accuracy can be increased up to 
machine precision without any additional programming 
effort. However, it would be misleading to infer that par- 
ticle mesh methods sacrifice accuracy in favor of speed, 
because due to the more advantageous scaling with parti- 
cle number (essentially N\ogN compared to iV 3 / 2 ) there 
will always be a critical number N*, such that the mesh 
method will be faster than the EwjaLd sum for particle 
numbers N > N*. See e.g. PetersenEZl for a discussion of 
the break-even value N* for PME. 



Error as a function of minimum image distance 

Instead of calculating the rms error for a complete con- 
figuration, it is also worthwhile to investigate it as a func- 
tion of the minimum image distance r between just two 
particles. This is a possibility to monitor the distance de- 
pendence of the accuracy for the various methods. Thus, 
we randomly created a pair of particles inside the sim- 
ulation box (again, L = 10) with given minimum image 
separation r and calculated the rms error from Eqn. (671). 
This was repeated for 5 x 10 4 separations equally spaced 
between and hy3L, which is the largest possible min- 
imum image separation. As this is done at constant a 
for each method, the real space contribution to the force 
always cancels when performing the difference in Eqn. 
(|37j), so this plot is only sensitive to the Fourier con- 
tribution and it is not necessary to specify a real space 
cutoff. A grid with Nm = 32 was chosen and the charge 
assignment order was set to P = 7. However, as can 
be seen from Fig. [|, different methods have their opti- 
mal operation point at different values of a. Therefore 
we found it more sensible to compare the different meth- 
ods at their individual optimal value of a, which can be 
obtained from Fig. 0. Although the curves in Fig. || cor- 
respond to a system which contains 100 particles (and 
not just 2), we believe that this has no influence on the 
optimal value of a, since e.g. the error formulas for the 
Ewald sum derived by Kolafa and Perramlia show, that 
the real space and the Fourier space contribution to AF 
display the same dependence on particle number. 

Note that the Coulomb problem in the given periodic 
geometry lacks spherical symmetry and due to the exis- 
tence of a grid also the translational symmetry is broken. 
So - strictly speaking - AF is not just a function of r but 
also depends on the orientation of the particles and their 
location within the box. This manifests itself in the fact 
that the measured points AF(r) do not collapse onto a 
single smooth curve but show some scatter. Since we are 
not interested in this effect, we averaged the scatter by 
binning 50 points together at one time and additionally 
performing a Gaussian smoothing (with width 0.1). This 
just makes the data easier to plot and digest. 

The result of this measurement is shown in Fig. [?]. 
Several interesting things can be observed: All algo- 
rithms produce their largest errors at small distances 
and get considerably more accurate at larger values of 
r - with one exception: The analytically differentiated 
SPME method almost immediately settles to a (compar- 
atively large) constant error. Since the only difference be- 
tween line 2 and 3 is the differentiation scheme, it must be 
the random force fluctuations discussed in the section on 
back- interpolation which are responsible for this effect. 
Note that PME at some distance gives better results than 
ik-diffcrentiated SPME. Also it is most surprising that 
at large distances PME and P 3 M give identical errors, 
although they differ considerably in the charge assign- 
ment scheme as well as in the employed Coulomb Green 
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r [C] 

FIG. 7. Distance dependency: The rms error AF as a func- 
tion of the minimum image separation r between two particles 
is shown for several mesh methods, which all share = 32 
and P — 7. Line 1 corresponds to PME, line 2 and 3 are an- 
alytically and ik-differentiated SPME respectively and line 4 
is ik-differentiated P 3 M. For each method a was individually 
set to its optimal value from Fig. ^: a.\ — 0.8, «2 = 0.89, 
a 3 = 0.92 and a 4 = 0.94. 

function. Finally, the ik-differentiated P 3 M method is 
most accurate for all distances. In this case this is not so 
much surprising, because the quantity Q from Eqn. (|29|), 
with respect to which P 3 M is optimized, is essentially the 
integral over any of these curves in Fig. [?], weighted with 
the probability density of the minimum image distance 
r. 



CONCLUSIONS 

Based on our theoretical considerations and the results 
of our numerical experiments we draw the following con- 
clusions: 

• The error of all Ewald calculations - be it the stan- 
dard Ewald sum or any particle mesh method - de- 
pends very sensitively on the Ewald parameter a. 
Hence, finding the optimal value of a is not merely 
an option but absolutely essential. Abstaining from 
a proper a-tuning results (at best) in wasting ac- 
curacy and (at worst) in the calculation of wrong 
forces or energies. For the standard Ewald sum 
and for PME there exist estimates for the real- and 
reciprocal space contribution to the error, which al- 
low an a priori determination of the optimal value 
for a, depending on the relevant system parame- 
ters. For the P 3 M method we will-tackle this prob- 
lem in a forthcoming publicationO. 

• At given computational effort the errors produced 
by different methods do not just vary marginally 



but by orders of magnitude, so putting some effort 
into this topic is certainly worth the trouble. 

• Generally, the total error is a combination of several 
contributions. If one of them dominates, there is no 
point in improving the other parts of the method. 
Assume for instance that we are using a finite dif- 
ference scheme for the derivative and that the over- 
all accuracy is actually limited by the discretization 
errors resulting from a low order charge assignment 
scheme. In this case, increasing the order J of the 
differentiation scheme would be useless. E.g., in 
Fig. U going beyond J = 7 would not yield any im- 
provement. If in this figure the assignment order P 
was 3 and not 7, it would even suffice to use J = 2 
(compare with Fig. ||). 

• Different methods scale differently with respect to 
their parameters. E.g., since the a-dependency of 
the error for PME is not the same as for SPME 
(there is not just a constant factor between them), 
the error curves can intersect (see Fig. The 
(almost trivial) consequence is that if one method 
is more accurate than another method in a specific 
region of the space of tuning parameters, this need 
not be the case in another region or even for all 
choices of parameters. 

• There exist many possibilities for combining the 
various parts of a mesh calculation - like charge 
assignment or differentiation scheme. This free- 
dom of choice can be exploited to suit ones particle 
mesh algorithm to already existing constraints in 
the complete simulation program. However, these 
combinations should always be tested thoroughly, 
since naive "improvements" can turn out to be dis- 
astrous (see line 2 in Fig. ||). There are, so to speak, 
several incompatible roads towards optimization, 
and one step away from a local optimum is in gen- 
eral a disimprovement. 

• If method A is at the same computational effort 10 
times more accurate than method B, this can be 
advantageous even if one is happy with the accu- 
racy of B: Almost surely method A - tuned down to 
the accuracy of B - will be faster than B, because, 
e.g., the number of mesh points could be reduced. 

• If one wants to use the continuum Green func- 
tion, a Lagrange interpolation scheme should be 
used. However, our tests show that a combination 
of the smooth spline interpolation with an appro- 
priately adjusted Green function - like in the P 3 M 
and SPME approach - should be preferred, since 
this can be made more accurate. We recommend 
the P 3 M approach, because it uses the analytically 
derived optimal influence function from Eqn. (|3l|), 
which minimizes the force errors, and - due to its 
smooth charge assignment - permits all investi- 
gated differentiation schemes. In particular, the 
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ik method is the most accurate implementation of 
the derivative (and comes most closely to the Ewald 
method, see Eqn. fll5|)), whereas the analytic differ- 
entiation (introduced by Essmann et. al.a originally 
for the SPME method) - although somewhat less 
accurate - is a fast and attractive alternative. 

• Compared to a standard Ewald sum, which uses 
the same number of k-vectors, all mesh algorithms 
are much less accurate (see Fig. ||). However, the 
accuracy which is actually needed in a simulation 
is typically not too large, since most simulations 
employ at the same time some kind of thermostat, 
and it is a waste of time to calculate the electro- 
static forces much more accurate than the random 
fluctuations of the thermostat. 

Note added in proof. After the submission of our 
paper T. Harden kindly brought to our attention a 
publication!!] where he performed a numerical compar- 
ison of the P 3 M to the SPME method, leading to results 
which are in agreement with our findings. 
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APPENDIX A: SOME COMMENTS ON UNITS 

Different people and communities prefer different con- 
ventions for units, especially if it comes to electrostatics. 
In this small appendix we present our choice. 

We write the Coulomb potential generated by a point 
particle with charge q located at the position ro as: 



r = 



<1 



(Al) 



Thus, its dimension is charge divided by length. In other 
words, if we measure all lengths in multiples of some unit 
length C and all charges in multiples of some unit charge 
C, the dimension of the electrostatic potential is C/C. As 
a consequence, the dimension of electrostatic energy is 
C 2 /£ and of electrostatic force is C 2 /£ 2 . 

In this article there is no need in specifying C or C, and 
if it comes to the final result (be it formulas or numbers), 
it can always be embellished with prefactors like I/Atteq. 

We give just one example: If one chooses C = A = 
lCT 10 m, C = e « 1.6022 x 10" 19 C and includes the 
standard-SI-prefactor I/Atteq, the numerical value of the 
original expression qtqz/r gives the force in units of 
2.3071 x 10" 8 N. 



Another common unit of force - especially among 
chemists - is kcal mol -1 A -1 . Obviously we have 



kcal 



4.186 x 10 3 J 



6.9510 x 10 _1U N 



mol A 6.022 1 7 x 10 23 lO" 10 m 

Thus, if one prefers to measure forces in units of kcal 
mol -1 A -1 , one only has to multiply the numerical value 
of the original qxqijr 1 by a factor of approximately 331.9. 



APPENDIX B: CONTINUUM GREEN 
FUNCTION AND LAGRANGE 
INTERPOLATION SCHEME 

In this appendix we show, how the implemented Green 
function and the charge assignment scheme are related to 
each other. More specifically, we demonstrate that the 
use of the continuum version of the Coulomb Green func- 
tion, as it appears in the conventional Ewald sum, sug- 
gests a so called Lagrange interpolation scheme, because 
this leads to a nice cancellation of certain discretization 
errors. We closely follow the notation of Hockney and 
EastwoodQ. 

We consider only the one-dimensional case. The elec- 
trostatic potential at position x' due to a unit charge 
residing at position x is not just a function of \x' — x\ 
but also depends on the distances of this charge from 
its neighboring mesh points. This artifact of the mesh 
can be quantified as follows: Let g(x) be the continuum 
Coulomb Green function and W p (x) = W(x — x p ) the 
charge assigned to mesh point p at position x p due to a 
unit charge at position x. The electrostatic potential at 
position x' can then be written as 



p 

E 

P =i 



4>{x')= Y,W p (x) g(x' - x p ) 



(Bl) 



where the sum is taken over all P mesh points to which 
the particle at position x contributed some fraction of 
its charge, i.e. the P mesh points which are closest to x. 
Taylor expanding g(x! — x p ) about (x' — x) gives: 



<f>(x') 



p 

E 

p=i 



w, 



pW E ^—^-g {n) (x - x') 



n=0 



(B2) 



It is possible to cancel the artificial terms in the n-sum 
(i.e. the ones which depend on x — x p ) up to order P by 
choosing the charge fractions W p (x) such that 

p 

J2w p {x) (x - xp)"- 1 = 5 hn , n=l,...,P (B3) 
P =i 

By induction with respect to n one can show that this 
may equivalently be expressed as 

p 

^Wp^x;- 1 =z"- 1 , n = l,...,P (B4) 

p=i 
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This system of P linear equations has a unique solution 
for the W p (x) since the coefficient matrix a;™ -1 is a Van- 
dermonde matrix for the distinct points Xi, . . . ,xp and 
hence has full rank. The W p {x) are thus polynomials of 
degree P—l. Since in particular Eqn. must be true 
at the mesh points, it follows 



X>p( 

P =i 



x q )x n ~ x =x n - x , n,q=l,...,P (B5) 



which - again due to the invertibility of x™ 1 - can only 
be true if 



W p (x q ) = 5. 



PQ 



p,q=l,...,P 



(B6) 



suffices to determine the polynomials 
referred to as 



Equation 

W p {x). They are referred to as the fundamerda, 
polynomials for the Lagrange interpolation problemEQ. 
PetersenO tabulates them for P = 3, . . . , 7 and their im- 
plementation is explained in detail. If one needs these 
assignment functions for higher values of P, one has to 
solve the syste m of linear equations (B4) or the interpo- 
lation problem (B6). 



APPENDIX C: SYSTEMATIC DIFFERENCE 
APPROXIMATIONS TO THE DIFFERENTIAL 
OPERATOR 

In this appendix we show, how mesh approximations 
for the differential operator d/dx can systematically be 
written as convex combinations of difference operators. 
In this way one can implement optimal combinations of 
these operators into the program- right-away, so an em- 
pirical tuning of the coefficientsLfEJ is no longer necessary. 
We describe the idea only for one dimension, the gener- 
alization to higher dimensions can be done easily via the 
Cartesian components. 

First we define the j th -neighbor centered difference op- 
erator Aj by 



(Aj f)(x) 



f(x + jh)-f(x-jh) 
2jh 



(CI) 



where h is the mesh spacing and x some mesh point. 
Applying this operator on a function / can be written as 
the convolution Dj ★ /, where Dj(x) is defined as 



DAx) 



8(x + jh) — 5{x — jh) 
2jh 



(C2) 



From the convolution theorem it follows that in Fourier 
space the derivative is given by Dj(k)f(k), where it is 
easily verified that the Fourier transform of Dj is 



.sm(jkh) 



(C3) 



(Note that in the limit h j. this reduces to the Fourier 
representation of d/dx, namely ik.) 



order J 


Cl 


C2 


C3 


1 


1 






2 


4/3 


-1/3 




3 


3/2 


-3/5 


1/10 


4 


8/5 


-4/5 


8/35 


5 


5/3 


-20/21 


5/14 


6 


12/7 


-15/14 


10/21 



<"4 



CO 



-1/35 
-5/63 



1/126 



TABLE I. Optimal form for the weighting coefficients of 



the J -order difference operator A'" 7 ' from Eqn. (C4) for 



several values of J. 



Since one can expect to achieve better approximations 
for the differential operator by using linear combinations 
of the difference operators Aj , we define a J th -order dif- 
ference operator by 



J 

£■ 



(C4) 



Using the Fourier representation of the differential and 
the j th -neighbor centered difference operator from Eqn. 
(C3), we demand 



J 

E 

j 



. sin(jkh) 



cji — jj— = ik + 0{{kh) 2J+1 ) (C5) 



or C J cos(jkh) = 1 + 0{{kh) 2J ) (C6) 

where the second equation follows from differentiating 
the first. Taylor expanding the cosine in Eqn. (C6) and 
equating coefficients gives J linear equations for the J 
unknowns Cj . The first few are given in Table (||) . 

Note that in the case of the 2 nd -order approximation 



the wei£ 
optimal 



iting (g, — g), which empirically was found to be 
, is reproduced. 



APPENDIX D: THE MODEL SYSTEM 

The rms error in the force for a system of 100 parti- 
cles randomly distributed in the simulation box is some- 
what sensitive to details of the generated configuration, 
e.g. the actual minimum distance. In order to make our 
measurements fully reproducible we decided to present 
our configuration as well. 

We found it easier not to list the particle positions 
but to describe the procedure which was used to gen- 
erate them. The coordinates of the 100 particles were 
constructed by first drawing 300 random numbers lZ n 
between and 1. If L is the box length then particle 
1 gets the coordinates (LIZi, LIZ2, LIZ3), particle 2 gets 
(LIZ4, L1Z 5 , L7Z 6 ) and so on. Moreover, particles with 
an even/odd number will get a positive/negative unit 
charge. 
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The choice of the random number generator is the fol- 
lowing: If a n is a positive integer, define its successor 
a n+ i via: 



- particle 



a n+1 := (1103515245 a» + 12345) mod 2 32 
Now define the pseudo random number lZ n G [0; 1[ by 
(a„ 65536) mod 32768 



(Dl) 



Tl n := 



32769 



(D2) 



Where "-i-" should denote an integer division which dis- 
cards any division rest. Choosing ao = 1 we obtained 
the sequence of random numbers 16838, 5758, 10113, . . ., 
of which the first 300 were used for positioning the parti- 
cles, e.g. - with L = 10 - the first particle has coordinates 
(5.138 . . . , 1.757 . . . , 3.086 . . .) and a negative unit charge. 
The smallest minimum image distance is approximately 
0.370264 and occurs between particle 46 and particle 98. 

Incidentally, we did not choose this random number 
generator because it is particulary good (it is not), but 
it is very easy to implement. Many C librarie s provide a 
function rand, which relies on Eqns. (Dl D2|). 



APPENDIX E: CHARGE ASSIGNMENT WITH 
SPLINES 

In this appendix we describe in a little more detail the 
procedure of charge assignment and present the charge 
fractions which are needed for a i2. th order assignment 



scheme a la Hockney and EastwoodQ (see Eqns. ( 18 . 1S|) ) ■ 
Let the units be chosen such that the grid spacing 
is 1. For any P consecutive mesh points there exists 
an interval X of length 1 such that the charge of a 
particle with coordinate x S T is distributed between 
these mesh points. By simple shifting we can assume 
this interval to be [— §,+g]. Then the P mesh points 
will lie at - E f 1 , + 1, . . ■, E f L - For P = 3 this 
is schematically shown in Fig. @. The charge fraction 
W p (x), which will be assigned to the mesh point x p , 
is related to the charge assignment function W(x) via 
W p (x) = W(x-x p ). 

The basic steps which have to be done for a particle 
with coordinate x (generally not in [— during a 

P th order charge assignment are thus: 



1. Define x to be the coordinate of the particle's near- 
est mesh point (if P is odd) or the midpoint be- 
tween the two nearest mesh points (if P is even). 

2. Find the P mesh points x p which are closest to x. 
They will be indexed by their relative position to 
x,sope {-£=1,-^1 + !,...,£=!} 



3. The fraction of charge which is assigned to each of 
these mesh points is given by W p (x — x). 





x 


T 


X 


x-axis 






X~» 

X 







-1 1 

FIG. 8. Schematic picture for a three-point charge assign- 
ment. The crosses are the mesh points and the lines indicate 
the (Wigner-Seitz) cell boundaries of each point (the mesh 
spacing is h = 1). All particles with x 6 [— §,+5] distribute 
their charge between the mesh points at -1, and +1 and the 
corresponding charge fractions, are W p (x), p G { — 1,0, +1}. 
(After Hockney and Eastwood!!) 



In this way the charge fractions are written as a func- 
tion of the separation x — x 6 [— 5, +5]. Hockney and 
Eastwood refer to the cases P = 1, 2 and 3 as NGP (near- 
est grid point), CIC (cloud in cell) and TSC (triangular 

shaped cloud) respectively. Generally, for P € {1, . . . , 7} 

(p) 

the charge fractions Wp (x) are given by the following 
polynomials: 
P= 1: 



W^(x) = l 



P = 2: 



W%{x) = -{l-2x) 



1 

2 l 

W%{x) = \{l + 2x) 



P = 3: 



W^{x) = -(1-4.t + 4x 2 ) 
8 

VF (3) (*) = i(3-4.x 2 ) 



W^{x) = ^(1 + 4.t + 4x 2 ) 
8 



1 



P = 4: 



W%(x) 


= 4^- 


6x + 


I2x 2 - 


8x 3 ) 


W%{x) 


= 4V 23 


-30x 


- 12x 2 


+ 2Ax 3 ) 


W%(x) 


= 4> 


f 30a; 


- 12x 2 


- 24x 3 ) 




= i( 1 + 


6x + 


I2x 2 + 


8x 3 ) 



P = 5: 



W (5 }{x) = — (1 - 8x+ 2Ax 2 - 32x 3 + 16z 4 
1 ' 384 v 



Wf\x) 



W K _{(x) = —(19 - 44a; + 24a; 2 + 16a; 3 - 16a; 4 ) 



W${x) 



1 

192 
1 



(115- 120a; 2 + 48a; 4 ) 

2 1CJ 



W${{x) = —(19 + 44a; + 24a; 2 - 16a; 3 - 16a; 4 ) 



384 



(1 + 8a; + 24a; 2 + 32a; 3 + 16x 4 
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P = 6: 



P = 7: 



W (6) , 

VV -5/2 



.a; 



WW 



-3/2 V* 
-1/2 V*- 



f r -1-1 / 



Vl/2( :E - 
"+3/2V . 



w (6) 

"-WW 



+5/2 V- 



a; 



3840 
1 

3840 
1 

1920 
1 

1920 
1 

3840 
1 

3840 



(1 - 10a; + 40a; 2 - 80a; 3 + 80a; 4 - 32a; 5 ) 
(237 - 750a; + 840a; 2 - 240a; 3 - 240a; 4 + 160a; 5 ) 
(841 - 770a; - 440a; 2 + 560a; 3 + 80a; 4 - 160a; 5 ) 
(841 + 770a; - 440a; 2 - 560a; 3 + 80a; 4 + 160a; 5 ) 
(237 + 750a; + 840a; 2 + 240a; 3 - 240a; 4 - 160a; 5 ) 
(1 + 10a; + 40a; 2 + 80a; 3 + 80a; 4 + 32a; 5 ) 



W^{x 

W^{x 
wP{x 

W$(x 



46080 
1 

23040 
1 

46080 
1 

11520 
1 

46080 
1 

23040 
1 

46080 



(1 - 12a; + 60a; 2 - 160a; 3 + 240a; 4 - 192a; 5 + 64a; 6 ) 
(361 - 1416a; + 2220a; 2 - 1600a; 3 + 240a; 4 + 384a; 5 - 192a; 6 ) 
(10543 - 17340a; + 4740a; 2 + 6880a; 3 - 4080a; 4 - 960a; 5 + 960a; 6 ) 
(5887 - 4620a; 2 + 1680a; 4 - 320a; 6 ) 

(10543 + 17340a; + 4740a; 2 - 6880a; 3 - 4080a; 4 + 960a; 5 + 960a; 6 ) 
(361 + 1416a; + 2220a; 2 + 1600a; 3 + 240a; 4 - 384a; 5 - 192a; 6 ) 
(1 + 12a; + 60a; 2 + 160a; 3 + 240a; 4 + 192a; 5 + 64a; 6 ) 
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